% This script is written to realistically simulate a tidal signals and
% resample it for night time only and try recover the amplitudes

%% Background spectra of slope -2
%define coefficients

clear

A = [1 -.999]; % slope  of -2
fs = 1/3600; % sampling every 1 minute

numsamples = (365 * 24 * 3600) * fs;

% create background red signal. Adjust the amplitudes to control the slope
% effect
background = filter(1,A,randn(numsamples,1));


%periods = [4.0000    4.8000    6.0000    8.0000   11.9672   12.0000   12.4210   12.6583   23.9345   24.0000   25.8910];

% The last period is to test aliasing effect to M2 signals from frequencies
% multiple of k* sampling frequency + M2 frequency. Here, the sampling
% frequency can be changed by the night sub sampling. In my experience, the
% aliasing effect from night sampling is not very high

Frq = [1/(18.6 * 365 * 24) 1/25.819 1/24.8421 1/12.4211 1/12.658]./3600;

%Frq = [1/24.8421 1/12.4211 1/11.9672]./3600;


% Add tidal lines .
Amp = ones([1,length(Frq)]);
noise_power = 0.1;


t = (1:numsamples)/fs;

y =     Amp(1)*cos(2*pi*t*Frq(1));

for i = 2:length(Frq),
    y = y + Amp(i)*cos(2*pi*t*Frq(i));
end;

y = y +   normrnd(0,noise_power,1,length(t));

sdata = y + background';

%% Recover the amplitudes
periods = 1./Frq;

% periods = [7.0000 17.0000 23.0000].*3600; some random periods

%make a model


m=[cos(2*pi*t(:)/(periods(1)) ) sin(2*pi*t(:)/(periods(1)) )];

for j=2:length(periods)
    m=[m cos(2*pi*t(:)/(periods(j)) ) sin(2*pi*t(:)/(periods(j)) )];
end

% now invert with the whole synthetic signals

[b, stats] = robustfit(m,sdata); % data with "red" noise
[a, stats] = robustfit(m,y); % data without "red" noise

a_amps = abs(complex(a(2:2:end),a(3:2:end)));
b_amps = abs(complex(b(2:2:end),b(3:2:end)));



fprintf('Expected White Red\n');
for i = 1:length(Frq),
    fprintf('%6.4f %6.4f %6.4f\n', Amp(i),a_amps(i),b_amps(i));
end;


%% Day - night separation

day_length = 3600 * 24 * fs;
night_length = round ( day_length * ( 6/24)); 


selected_blocks = false ( 1, length(y));

for i = 1:length(y) / day_length,
    
    selected_blocks (((i-1) * day_length) + 1 : ((i-1) * day_length) + night_length) = true;
    
end;

y_night = y(selected_blocks);
s_night = sdata(selected_blocks);
m_night = m(selected_blocks, :);
t_night = t(selected_blocks);


[c, cstats] = robustfit(m_night,y_night); % data without "day"
[s, sstats] = robustfit(m_night,s_night); % data without "day"


c_amps = abs(complex(c(2:2:end),c(3:2:end)));
s_amps = abs(complex(s(2:2:end),s(3:2:end)));


fprintf('Period Expected All Night RedN\n');

for i = 1:length(Frq),
    fprintf('%5.2f %6.4f %6.4f %6.4f %6.4f\n', 1./(Frq(i)* 3600),  ... 
        Amp(i),a_amps(i),c_amps(i), s_amps(i));
end;


%% plot the original and recovered data

% Get aplitude and phase 
P = angle(complex(c(3:2:end),c(2:2:end)));
A = abs(complex(c(2:2:end),c(3:2:end)));
if exist ('yy', 'var'),
    clear yy;
end;
for i = 1:length(P),
   yy(i,:) =   A(i) * sin(2*pi*t./periods(i) + P(i));
end;

plot(t, y);
hold on
plot(t_night, y_night, 'r.');
plot(t, sum(yy),'k');


